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We consider the problem of calculating the large number of Wick contractions necessary to com- 
pute states with the quantum numbers of many baryons in lattice QCD. We consider a constructive 
approach and a determinant-based approach and show that these methods allow the required con- 
tractions to be performed in computationally manageable amount of time for certain choices of 
interpolating operators. Examples of correlation functions computed using these techniques are 
shown for the quantum numbers of the light nuclei, 4 He, 8 Be, 12 C, 16 and 28 Si. 



The ab initio approach to nuclear physics from the underlying theory of the strong interactions, Quantum Chro- 
modynamics (QCD), is hampered by the many body nature of the nuclear problem. In principle, QCD and the 
electroweak interactions give rise to all the rich and complex phenomena of nuclear physics, yet it is only recently 
that the first QCD studies of multi-baryon systems have appeared [IHH] • The reason for this is twofold. Firstly, the 
Monte-Carlo evaluation of correlation functions of multi-baryon systems converges slowly, requiring a large number 
of measurements before the necessary precision is reached (this issue will not be addressed here). Secondly, systems 
with the quantum numbers of many nucleons and hyperons are complex many-body systems with complicated spectra 
and there are a multitude of physically relevant states that can be studied in QCD. Even for a given set of quantum 
numbers, additional complexity appears at the quark level; the number of Wick contractions required to construct 
systems for large atomic number grows factorially, scaling as n u !nd!n s ! where n ai j s are the numbers of up, down, 
and strange quarks required to construct the quantum numbers of the state in question. In many situations, this is 
a naive counting as there are many cancellations and contributions that are identical. However, the a priori identi- 
fication of these simplifications is a non-trivial task. In addition to the problem of Wick contractions, the number 
of terms in the interpolating fields of multi-nucleon systems also typically grows exponentially with the size of the 
system. This potentially more serious problem is similar in nature to the problem of the exponential growth of nuclear 
wave-functions faced in nuclear structure calculations where phenomenological potential models describing the low 
energy nucleon-nucleon interactions are used. 

In this paper, we present a systematic method for the construction of nuclear interpolating fields for multi-baryon 
systems in lattice QCD (LQCD) (see Ref. [5] for related work). We demonstrate that the Grassmannian nature of the 
quark fields can be used to our advantage, in some cases resulting in particularly simple nuclear interpolating fields. 
In addition, we present two approaches that ameliorate the cost of contractions, the most efficient of which scales 
only polynomially in the number of quarks involved in the contraction. Using these methods we compute LQCD 
correlation functions with the quantum numbers of the light nuclei, 4 He, 8 Be, 12 C, 16 and 28 Si, demonstrating that 
correlation functions relevant to the study of nuclei in QCD can be constructed. 



In order to calculate nuclear correlation functions, we first need to construct quark level nuclear interpolating fields. 
This is, in principle, straightforward and, in practice, it resembles the construction of quark model wave-functions for 
baryons |10) . A general quark-level nuclear interpolating field with atomic number A containing n q — 3 A quarks has 
the form 



where the q ai are the quark fields, the are generic indices which combine the colour, spinor, flavour, and spatial 
indices of the quark and a is a compound index representing the n g -plet a\,a2 ■ ■ ■ a n . Given that calculations are 
performed on a discrete lattice, the spatial degrees of freedom are finite and countable, and as a result we can use 
an integer index to describe them. Here the quark fields are all at the same time t. The index h on the nuclear 
interpolating field is a set of quantum numbers that identify the nuclear state, including its momentum, angular 
momentum, isospin and strangeness. The Grassmannian nature of the quark field dictates that the tensor w^' 2 " q 
is totally antisymmetric under the exchange of any two indices. If the indices can have a total of N possible values, 
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then, ignoring the detailed flavour structure, the total number of non-vanishing terms in the above sum is 



(N-n q )\ 



(2) 



However, many of these terms correspond to permutations of the quark fields, and the total number of unique terms 
(terms that are not a permutation of any other term) is 

N\ 

(3) 



n q \{N-n q )\ ■ 

For a generic spatial structure of the interpolating field, N corresponds to the total number of quark degrees of freedom 
on a time-slice (N = 12L 3 , where L is the spatial dimension of the lattice) and one may be discouraged by the feasibility 
of the task of building quark level nuclear interpolating fields. However, a number of simplifying factors are omitted in 
the above discussion. As we consider interpolating fields with definite transformation properties under the symmetries 
of QCD, large numbers of terms in the nuclear interpolating field vanish. The first major reduction comes from the 
fact that only colour singlets need to be considered. In addition, considering only interpolating fields of definite parity, 
angular momentum 1 , isospin and strangeness, forces several elements of the tensor 2 "'to vanish. Finally, the 
most drastic reduction of the non-zero tensor elements can be achieved using simple spatial wave-functions. At this 
time, having recognised that only a small fraction of the terms in the sum of Eq.[l]are non-zero, as well as the fact that 

the tensor w h ' is totally anti-symmetric, we can introduce the reduced weights w h which are the 

minimal set of non-zero numbers required to completely describe the interpolating field. The n 9 -plet (ai, a 2 • • • a nq ), 
is an ordered list of indices that represents a class of terms in Eq. [I] that are all permutations of each other. The 
index k on the reduced weights enumerates the number of classes that the tensor w h ' 2 nq decomposes into. With 
these reduced weights, Eq. [T]can be re- written as 

x h = E w { : i ' n ^' k J2 iin ° ?(«,)?(«,) • • • «k ) < ^ 

k=l i 

where N w is the total number of reduced weights, l represents the ?ig-plet (21, £2 * * * ^n q 

) and e 41 '* 2 '"' ' tn i is a totally 

anti-symmetric tensor of rank n q with 

1,2,3,4,- ,?», _ 1 

The above expression is the simplest form of the quark-level nuclear interpolating field and is completely described 
by the reduced weights. As an example, using a single point spatial wave function, the numbers of terms contained 
in the simplest interpolating fields for the proton, deuteron, 3 He and 4 He, are N w = 9, 21, 9, and 1, respectively. 



A. Hadronic Interpolating Fields 

Having now written down a general nuclear interpolating field with quantum numbers h, we need to calculate the 

reduced weights w~ h 12 n * in an efficient manner. In principle, this can be achieved directly from quark fields 
by imposing the desired transformation properties. However, in certain cases, it is advantageous to proceed by first 
constructing hadronic interpolating fields from which the quark interpolating fields arc derived. 

The hadronic interpolating fields assume a form analogous to that of the quark interpolating fields. The baryons 
that make up the nucleus are also fermions, hence the general structure outlined above can be directly transcribed 
here. In terms of baryons, a nuclear interpolating field of a nucleus of atomic number A is 

M w 

M h = W^> b2 ~ bA) ^' i2 '"' M B{K)B{b i2 ) ■ ■ ■ B(b lA ) , (5) 

fe=l i 

where M w is the number of hadronic reduced weig are baryon interpolating fields and the bi are 

generic indices that includes parity, angular momentum, isospin, strangeness, and spatial indices. Unlike the quark 



For simplicity, we refer to the irreducible representation of the lattice symmetry group as angular momentum. 
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fields which are fundamental degrees of freedom, the baryon interpolating fields are composite objects, hence there is 
a large number of such interpolating fields for a given set of quantum numbers. For simplicity, as well as efficiency 
of the resulting nuclear interpolating fields, we will use a single interpolating field per baryon, selected so that it has 
good overlap with the single baryon ground state, as well as being comprised of a small number of quark level terms. 
The utility of the above form of the nuclear interpolating fields is twofold. Firstly, it allows us to derive the reduced 
weights we need for Eq. [4] Secondly, interpolating fields that are derived starting from Eq. [5] may have better overlap 
with the nuclear ground states as it is well-known that hadronic degrees of freedom provide a successful description 
of much of nuclear physics. 

The calculation of the reduced weights, Wv 6l ' &2 bA \ in the hadronic interpolating field is straightforward. It 
amounts to combining individual hadrons of given quantum numbers to build a multi-hadron state of definite parity, 
angular momentum, isospin, and strangeness. This construction can be readily automated and can be performed 
recursively using the known Clebsch-Gordan coefficients of SU(2) for both the spin and isospin (or SU(3) flavour if so 
desired). In principle, one can use all the octet and decuplet baryons in Eq.pl however, for most practical purposes, 
restricting to the positive parity octet baryons is sufficient. For example, for^l = 2, I = J = 0, S = —2, if we restrict 
the spatial wave-function to single point, there are three simple hadronic interpolating fields 



v3 

and, 

l - [s°M - s-V - s ^ + s~V] , 

where and represent the spin up and down polarisations of the baryon, B, respectively. In this example, the 
reduced weights can be directly read off from these equations. 

We have written a C++ symbolic manipulation program that generates the hadronic reduced weights using the above 
approach. In Ref. [5], we have used this to produce a complete basis of orthonormal interpolating fields with spatial 
wave-functions restricted to a single point for all nuclei up to A = 4 and have also constructed a selection of states for 
A > 4. Generically for larger A, more complicated spatial wave-functions are required because of the Pauli exclusion 
principle, resulting in an exponential growth of the number of possible interpolating fields as A increases (this reflects 
the problem faced in nuclear structure calculations as A becomes large). In certain cases, the Grassmannian nature 
of the quark fields is also advantageous, drastically reducing the number of non-zero reduced weights. Making use of 
this feature, we have been able to find particularly simple wave- functions for systems as large as A = 28. 



B. Quark Interpolating Fields 

The reduced weights of the quark interpolating fields of Eq. [3] can be calculated by equating the two forms of the 
nuclear interpolating fields 

M w 

M h = ^^ 1 '^"- 6 ^^ e i »^'"-^B(^JB(6 ia )... J B(6 iA ) 

k=l i 

N w 

= ^„.; ' • ■ >^g( ail )g(a i3 ) • ■ -q(a inq ), (6) 

k=l i 

and replacing the baryon objects by their quark interpolating fields. A single baryon interpolating field is written in 
terms of quark fields as 

k=l i 

where Nb(i>) is the number of terms in the single baryon B(b) interpolating field. For single baryon interpolating fields, 
the weights, 1 ]j( ai ' a2 ' a3 ' > ' k ; have been presented in [11] (the colour factors necessary for our formulation are not included 
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in Ref. [IT] but can be trivially added). The process of deriving the reduced weights w^ 1 ^ 2 a ™ 9 ' from Eq. [ej can 
be automated and we perform it within our symbolic manipulation program. An interesting feature that arises from 

the calculation of the reduced weights tyf" 1 '" 2 °" nq ^ l,k j s that if we restrict ourselves to simple spatial wave-functions 
making use of only few spatial points, then the expected exponential growth of the number of terms in the nuclear 
interpolating field is eliminated. A careful selection of the spatial wave-functions used can eliminate this problem, in 
principle for arbitrarily large nuclei. However, restriction to a small number of quark degrees of freedom also makes 
it impossible to construct certain states (an example is presented in Ref. |8j where the two baryon symmetric flavour 
octet was found to be inaccessible). 

III. TECHNIQUES FOR MULTI-BARYON CONTRACTIONS 

In this section, we consider how the interpolating fields constructed in the previous section can be used to generate 
the correlation functions of multi-baryon systems. A general multi-hadron two point function is given by 

(•A/? (t)tf£ (0)> = i J VUVqVq TV? (t)Af£ (0) e- 5 « CD , (8) 

where Sqcd and Z are the QCD action and partition function respectively, and T>U, T>qDq are the gluon and quark 
field integration measures respectively. We have also introduced explicit dependence of the interpolating fields on the 
Euclidean time separation, t, and consider a two point function with different creation and annihilation interpolating 
fields with commensurate quantum numbers. For a given choice of the interpolating fields, it is straightforward to 
perform the Grassmann integral over the quark fields and re-write the correlation function in terms of the quark 
propagators. However, for an efficient calculation of the two point function we need to be mindful of the structure of 
the interpolating fields. 

One successful class of interpolating fields for two or more hadron systems is one in which a plane wave basis at the 
level of the hadronic interpolating fields is used. This amounts to projecting the individual hadrons comprising the 
multi-body system to definite momentum states, while preserving the spatial transformation properties of the overall 
multi-hadron system [TJ H^MTS] . In this case, the complexity of the spatial wave-function is such that the number of 
terms contributing to Eq. [4] is rather large and hadronic interpolating fields have to be used in order to build the 
desired two point function. Constructing these types of interpolating fields both at the source and the sink becomes 
computationally expensive because a large number of quark propagators that are required. Nevertheless, this method 
has been employed for meson-meson and multi- meson spectroscopy |18H21j . For the case of multi- meson systems, 
special contraction methods were required [HI [201 H2] • For multi-baryon systems, the problem is more complex and 
will be the subject of further investigations. A further approach is to consider correlation functions in which the quark 
creation interpolating fields (source) have simple spatial wave-functions with few degrees of freedom (for example, 
restricted to a few spatial locations), while using a plane wave basis for the hadronic interpolating fields at the sink. 
Finally, as we shall discuss below, sufficiently simple nuclear interpolating fields exist, where the number of terms 
contributing in Eq. [4] is small and factorization into hadrons is not computationally necessary. 

A. Hadronic blocks 

The quark propagator from a single source point, xq — (xo,0), can be used to construct baryon building blocks 
with quantum numbers b and momentum p, as: 

N B(b) 

B a b im (p, t; x ) = J2 eip " E *i Cl ' C2 ' C3) ' fc J2 ^' l2 ' l3 S(cn , x- a 1 ,x )S(c h , x; a 2 , x )S(c l3 ,x; a 3 , x ) , (9) 

x k=l i 

where S(c, x, i; a, xq, 0) is the quark propagator from xo to x — (x, t) and ai are the remaining combined spin- 
colour-flavour indices. In this notation, the sink indices are kept to the left of the source indices and the spatial 
indices are displayed explicitly as they play an essential role in the construction of the block. This baryon block 
corresponds to the propagation of an arbitrary three-quark state from the source to the sink where it is annihilated 
by the prescribed baryon interpolating field. As discussed above, we have chosen to momentum project these blocks 
at the sink to a given momentum p to allow control of the total momentum of multi-hadron systems, although this 
is not necessary and other forms of blocks can be envisaged. 
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We can generalise these blocks to allow the quark propagators to originate from different source locations, 

(1) (2) 

Xq , Xq , . . ., as necessary, using 

N B(b) 

Sr Q2 ' Q3 (p, t; s 1 ,s 2 ,s 3 ) = J2 e 4p ' X E *i Cl ' C2 ' C3) ' fc E ^'^(c^x; a u x^)S(c i2 , x; a 2 , a^JSfo , x; a 3 , 4* 3) ) , 

x k=l i 

(10) 

where the x label the source locations. These blocks can be further generalised to allow for non-trivial single hadron 
spatial wave-function at the sink, but we will not consider this case further. It may also be advantageous to consider 
more complicated multi-hadron blocks similar to those implemented in Ref. [2] although the storage requirements 
grow rapidly with number of baryons in the block. 

B. Quark- hadron contractions 

Using the building blocks described above, we can consider correlation functions in which quark level interpolating 
fields are used at the source and their hadronic counterparts are used at the sink. The contractions are performed by 
iterating over all combinations of source and sink interpolating field terms and connecting the source and sink with 
the appropriate sets of quark propagators. For a given pair of source and sink interpolating field terms, this amounts 
to selecting the components dictated by the source quark interpolating field from the product of blocks dictated by 
the hadronic sink interpolating field. The Wick contractions are implemented by performing this selection in all 
possible ways. This proceeds by taking the first hadron in the hadronic wave-function at the sink, replacing it by the 
appropriate hadron block and selecting the three free indices in all possible ways from the pool of indices dictated 
by the source quark interpolating field, keeping track of the appropriate permutation sign. Following this, the second 
baryon component in the hadronic (sink) interpolating field term is replaced with the appropriate block and the free 
indices are contracted with the remaining free indices in the source quark interpolating field term in all possible ways. 
These first steps are illustrated in Fig. [I] and the procedure continues until all hadrons in the sink interpolating field 
term have been contracted, necessarily using all available quark indices at the source. The result is then multiplied by 




(a) (b) 

FIG. 1: Illustration of steps one and two of the quark-hadron contraction method. The small circles in the left hand of 
the figures correspond to the quarks in the source interpolating field while the large squares and lines extending from them 
correspond to the hadronic blocks. 

the weights of the source and sink terms under consideration and added to the correlation function. The contraction 
is complete after all combinations of source and sink interpolating field terms have been considered. The process 
described here is independent of the the source and sink interpolating fields and can be applied to any correlation 
function. Further reductions of the total cost of the algorithm may be possible by studying the symmetry properties 
of a particular pair of source-sink interpolating fields. However, such reductions are not generic, hence we do not 



() 



consider them further. The procedure described has been used to perform the contractions needed for the large class 
of interpolating fields considered in the study of the spectrum of hyper- nuclei up to A — 5 in Ref. [U [23] . 

For large numbers of baryons (A > 8 for protons and neutrons alone) , it is necessary to use multiple source locations 
because of the Pauli exclusion principle. In this case, the generalised blocks in Eq. [lUJcan be used with the algorithm 
presented above. 



C. Scaling 

From the above description, it is clear that this algorithm will in general scale as 



v/, ' v '"'frF ! 



where A is the atomic number and M w and N w are the number of terms in the sink and source interpolating fields 
respectively. In addition, the fact that the hadron blocks are completely anti-symmetric under all quark exchanges 
has been taken into account. If we also take into account that the strong interactions are flavour-blind and consider 
only octet baryon building blocks, this reduces to 

M W -N W 2 H ^_i , (12) 

where n s o and n\ are the number of E° and A baryons in the hadronic interpolating field and the factor in the 
denominator arises because all octet baryons have two quarks of the same flavour except from the S° and A. This 
algorithm can be efficiently implemented and is computationally feasible for small systems, A < 10. As an example 
of this method, a 4 He two point correlation function can be computed in ~ 0.8 seconds per time slice on a single core 
of a Dual Core AMD Opteron 285 processor. 



IV. MULTI-BARYON CONTRACTIONS WITH DETERMINANTS 



For larger atomic number, A > 10, alternative methods are required to perform the contractions in a computationally 
feasible manner. It is straightforward to see how this can be done by examining the two point functions above and 
making use of Wick's theorem [24] . The numerator of Eq. [8] before the integration over the gauge fields is performed 
is given by 

K N w , , , , N , 

_/(o 1 ,O a — O, % ),k „(ai,a 2 ---a„ ),fc 



J fc'=ifc=i 
^2^ Ja, ''' J '" ,eil,<a '"' ,< " 4 9(°i».)-'-«( a i a )9( a iJ X q(a ll )q(a t2 )---q(a inq ), (13) 

j i 

where the primed and unprimed indices are associated with the sink and source interpolating fields, respectively and 
are composite colour, spinor, flavour and spatial indices and [..-]u indicates the value of the enclosed expression on 
a fixed gauge field. The Grassmann integral over quark fields can now be performed, resulting in the replacement of 
the qq pairs by elements of the quark propagator. 

K N w , , , 

k' = l k=l 

£ e^>"' S(a' h ; a n )S(a' 32 ;a l2 ) ■ ■ ■ S(o£ ; a lnq ) , (14) 



where S e ff[U] denotes the pure gauge part of the QCD action together with the logarithm of the determinant of the 
Dirac matrix. The above expression of Wick's theorem, can be written in terms of the determinant of a matrix G 
whose matrix elements are given by 

cm- ^ -I s{ - a p a ^ for a 'j e a ' and a * e a n 5^ 

^ a ' a ^-\5 0> otherwise ' [Lb) 



FIG. 2: Illustration of the quark determinant level contraction with each of the three sub-blocks listing the up, down and 
strange quarks of the wave-function term respectively . Within each flavour of quark, all possible contractions are performed 
by forming a determinant of the matrix of quark propagators as described in the text. 



where, as before, a' = (a' 1; a' 2 ■ ■ ■ o! n ) and a = (a±, a% • • • a nq ). Note also that the non-trivial block of the matrix 
G(a';a) is of size n q x n q , hence for computing its determinant we only need to consider this block. For this reason, 
in the following discussion, the matrix G(a'; a) denotes only this small non-trivial block. 
Making use of this definition, the full nuclear correlation function can be written as 

N' w N m 

(^(^(0)>4 We " S< " ^^< (a ' 1 ' aW - ) ^'4 Q - Q2 '-- a - ) ' fe xdetG(a';a). (16) 

The determinant of a matrix of size n q can be evaluated in n q operations (for example via LU decomposition) instead 
of the naive n q \ operations, so making use of this representation of the nuclear correlation function is numerically 
advantageous. Furthermore, because of the flavour-blindness of the strong interaction, the matrix G(a';a) is block 
diagonal, as a result the determinant calculation breaks into a product of smaller determinants, one for each flavour. 



Given the reduced weights determined above and appropriate quark propagators, the implementation of Eq. 16 is 
very fast, scaling polynomially with the number of terms in the source and sink quark level interpolating fields as well 
as the number of quarks per flavour. The total cost of this form of contractions scales as 

nlnlnl x N' W N W (17) 

where N' w N w are the number of terms in the sink and source quark interpolating fields respectively. 2 As a result, if 
we can construct interpolating fields with sufficiently small number of terms, correlation functions with a very large 
atomic number A can be constructed. 



V. NUCLEAR CORRELATION FUNCTIONS 



We have performed preliminary studies to investigate the numerical efficiency of these methods. Results for the 
quark-hadron approach have been presented in Ref. [5] and here we focus on the determinant-based approach. Cal- 
culations are performed on an ensemble of gauge configurations generated with a tadpole-improved Luscher-Weisz 



2 The expectation of polynomial scaling of contractions was noted by D. B. Kaplan in Ref. 25 . However, the scaling of N w and N' u 
with the atomic number A is generically exponential. 



gauge action and a clover fermion action with tadpole-improved tree-level clover coefficient. The gauge links entering 
the fermion action are stout smeared, with p = 0.125. Three flavours of quarks with masses corresponding to the 
physical strange quark mass were used. The lattice spacing, a ~ 0.145 fm, and the dimensions of the lattice are 
L 3 x T = 32 3 x 48 corresponding to a physical volume of (4.6 fm) 3 x 7.0 fm (further details will be presented elsewhere 
[2"6"]h We have performed a large number of measurements from spatially distinct sources on an ensemble of a 0(250) 
gauge configurations well separated in HMC evolution time. All calculations are performed in double precision and 
care is taken to preserve the dynamic range of correlation functions by rescaling quark propagators before contractions 
are performed. 

In Fig. [3j the logarithms of correlation functions are shown for correlators with the quantum numbers of the light 
nuclei, 4 He, 8 Be, 12 C, 16 and 28 Si. Error bars that reach the lower axis of the plots indicate that the correlator has 
fluctuations that are negative at one standard deviation. The extracted energies for each of the atomic number A < 20 
systems are consistent with a system of A nucleons but with large uncertainties at present (for 28 Si, no flattening of 
the effective mass is seen before the signal is lost) . Given the large number of near-threshold energy levels expected in 
these complex nuclear systems (see Ref. [8] for an example for A = 4) , a clean extraction of the ground state binding 
energies of these systems is beyond the current work. In addition, the baryon number density of the larger systems 
(0.3 fnU 3 for 28 Si) is substantial and volume effects are expected to be significant. It will be necessary to use larger 
volumes, increase greatly the statistical precision and improve the interpolating operators that we have used in order 
to obtain the binding energies and excitation spectra of these systems. Nevertheless, this study demonstrates the 
computational feasibility of lattice QCD calculations of light nuclei. 

VI. CONCLUSIONS AND OUTLOOK 

In this work, have discussed a systematic way of constructing interpolating fields for multi-baryon systems. In 
addition, we have investigated the issue of performing the Wick contractions with these interpolating fields in lattice 
QCD (see Ref. [9] for related work). We have shown that there are approaches that enable calculations of systems 
with very large number of nuclei that are computationally feasible and demonstrated their effectiveness by calculating 
correlators with baryon number up to A = 28. Given the expected spectra of such complex systems, significant 
advances are required in order to extract ground state energies from these correlators. Finally, the methods described 
here may prove useful in calculations of QCD at non-zero baryon density, where projection on to a given baryon 
number is required. 
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FIG. 3: Correlation functions for nuclear systems, 4 He, 8 Be, 12 C, 16 and 28 Si. In each row, correlators based on both 
smeared-point (SP) and smeared-smeared (SS) quark propagators are shown. 



